
********************************************************************************
*** Do file to replicate paper results 
*** Paper Title: Farm-Level Agricultural Productivity and Adaptation to Extreme Heat
*** Authors: Joaquin Mayorga, Alexis H. Villacis, Ashok K. Mishra 
********************************************************************************

* Clear environment 
clear all 

* Define paths 
global root_dir "C:\Users\jmayorga\Dropbox\AJAE Submission"
global data "${root_dir}\Data"
global shape_files "${root_dir}\Shape files"
global output "C:\Users\jmayorga\Dropbox\Apps\Overleaf\AJAE Accepted Tables and Figures"


* Load farm-level dataset 
use "${data}/final_dataset.dta", clear 

* Drop observations 
keep if d_final==1 

********************************************************************************
*** Define globals of variables  
********************************************************************************

* Inputs 
global inputs_ln "ln_ha_planted ln_labor_family ln_labor_hired ln_inorgfert_kg ln_herb_kg ln_pest_kg"
global inputs_dummies "labor_family_zero labor_hired_zero inorgfert_kg_zero herb_kg_zero pest_kg_zero"
global inputs "ha_planted labor_family labor_hired inorgfert_kg urea_kg npk_kg herb_kg pest_kg"


* Globals 
global soilall "i.sq1 i.sq2 i.sq3 i.sq4 i.sq5 i.sq6 i.sq7"
global temp "temp_avg temp_avg2"
global tempbins_linear "sh525less sh530 sh535plus"
global ltempbins_linear "lsh525less lsh530 lsh535plus"
global precip "rf rf2 rf_lag rf2_lag"
global rf "precip_daily_0308 precip_daily_0308_sqrd"
global lrf "rf rf2"
global lrf "rf_lag rf2_lag"
global rfz "rf_zscore rf2_zscore"
global hhchar "age age2 edu female lnhhsize_num"
global soilhalf "sq1moderateplus sq2moderateplus sq3moderateplus sq4moderateplus sq5moderateplus sq6moderateplus sq7moderateplus"
global soil "sq1 sq2 sq3 sq4 sq5 sq6 sq7"
global farmer "age age2 female edu"
global irrig "shareirrigated"

********************************************************************************
*** Figure 1: Location of EAs in sample 
********************************************************************************

preserve 
	keep lat_dd_mod lon_dd_mod ea year
	collapse (first) lat_dd_mod lon_dd_mod, by(ea)
	tempfile coordinates 
	save coordinates, replace
	shp2dta using "${shape_files}/nga_admbnda_adm1_osgof/nga_admbnda_adm1_osgof_20161215", database(nigeria_states_map) coordinates(coord) genid(id_ent) replace
	use nigeria_states_map, clear
	gen south = 0 
	replace south = 1 if admin1Name=="Abia"
	replace south = 1 if admin1Name=="Anambra"
	replace south = 1 if admin1Name=="Ebonyi"
	replace south = 1 if admin1Name=="Enugu"
	replace south = 1 if admin1Name=="Imo"
	replace south = 1 if admin1Name=="Akwa Ibom"
	replace south = 1 if admin1Name=="Bayelsa"
	replace south = 1 if admin1Name=="Cross River"
	replace south = 1 if admin1Name=="Delta"
	replace south = 1 if admin1Name=="Edo"
	replace south = 1 if admin1Name=="Rivers"
	replace south = 1 if admin1Name=="Ekiti"
	replace south = 1 if admin1Name=="Lagos"
	replace south = 1 if admin1Name=="Ogun"
	replace south = 1 if admin1Name=="Ondo"
	replace south = 1 if admin1Name=="Osun"
	replace south = 1 if admin1Name=="Oyo"
	label define south 0 "North" 1 "South"
	label values south south 
	spmap south using coord, id(id_ent) point(data(coordinates) y(lat_dd_mod) x(lon_dd_mod) fcolor(green) size(small)) ///
	clmethod(custom) clbreaks(-1 0 1) legcount fcolor(gray*.2 gray*.6) ///
	legend(ring(0) pos(7) col(1) size(medium) order(2 "North" 3 "South"))
	graph export "${output}\Figures\fig1_ea_locations.pdf", replace
	graph close 
restore 

********************************************************************************
*** Table: Descriptive statistics (need to update paper table) 
********************************************************************************

global hhchars "age edu female hhsize_num"
global prod "output_2010usd output_usdperha"
global inputs "ha_planted labor_hired labor_hired_ext fertilizer fertilizer_cropareapc urea_kg npk_kg herbicide herbicide_cropareapc herb_kg pesticide pesticide_cropareapc pest_kg"
global dstats "$hhchars $prod $inputs"

estpost summ $dstats 
est store all

esttab all,  ///
  cells( "mean(fmt(2) label(Mean)) sd(fmt(2) label(SD)) min(fmt(2) label(Min.)) max(fmt(2) label(Max.))" ) label ///
  stats(N, fmt(%12.0fc))
eststo clear

esttab all using "${output}/Tables/tab1_descriptive_stats.tex", replace  ///
  cells( "mean(fmt(2) label(Mean)) sd(fmt(2) label(SD)) min(fmt(2) label(Min.)) max(fmt(2) label(Max.))" ) ///
  label booktabs collabels(none) gaps ///
  stats(N, fmt(%12.0fc) labels("\midrule Observations")) ///
  title("Household and Farm Descriptive Statistics")
eststo clear

********************************************************************************
*** Figure: Percentage of planted area by month (still need to check sample)
********************************************************************************

preserve
	forvalues j=1/12{
		summ year if !missing(cropareasqm`j'_pc), d 
		rename cropareasqm`j'_pc cropareasqm_pc_`j'
	}

	collapse (mean) cropareasqm_pc_1-cropareasqm_pc_12, by(year)
	reshape long cropareasqm_pc_, i(year) j(month)
	rename cropareasqm_pc_ share_planted
	replace share_planted = share_planted*100
	collapse (mean) share_planted, by(month)
	label var share_planted ""
	label define month 1 "Jan" 2 "Feb" 3 "Mar" 4 "Apr" 5 "May" 6 "Jun" 7 "Jul" 8 "Aug" 9 "Sep" 10 "Oct" 11 "Nov" 12 "Dec"
	label values month month 
	graph bar share_planted, over(month) bar(1, color(green) lcolor(black)) ytitle("Percentage of annual area planted") 
	graph export "${output}\Figures\fig3_share_planted_by_month.pdf", replace
restore

********************************************************************************
*** Table: Weather statistics during the growing season and historical comparison (need to update paper table)
********************************************************************************

global tempbins_current 	"sh525less sh525 sh530 sh535plus"
global tempbins_long_term 	"sh525less_mean sh525_mean sh530_mean sh535plus_mean"
global tempbins_z_scores 	"sh525less_zscore sh525_zscore sh530_zscore sh535plus_zscore"
global climate_dstats 		"$tempbins_current rf $tempbins_long_term rf_mean $tempbins_z_scores rf_zscore"  

estpost summ $climate_dstats 
est store climate_dstats 

esttab climate_dstats ,  ///
  cells( "mean(fmt(2) label(Mean)) sd(fmt(2) label(SD)) min(fmt(2) label(Min.)) max(fmt(2) label(Max.))" ) label ///
   refcat(sh525less "Panel A: Weather during the growing season (2010, 2012, 2015)" ///
		sh525less_mean "Panel B: Weather over previous 25 years (1984-2009)" ///
		sh525less_zscore "Panel C: Z-scores (2010, 2012, 2015)" , nolabel) ///
///
  stats(N, fmt(%12.0fc))
eststo clear

esttab climate_dstats using "${output}/Tables/tab2_weather_statistics.tex" , replace ///
  cells( "mean(fmt(2) label(Mean)) sd(fmt(2) label(SD)) min(fmt(2) label(Min.)) max(fmt(2) label(Max.))" ) label ///
   refcat(sh525less "Panel A: Weather during the growing season (2010, 2012, 2015)" ///
		sh525less_mean "Panel B: Weather over previous 25 years (1984-2009)" ///
		sh525less_zscore "Panel C: Z-scores (2010, 2012, 2015)" , nolabel) ///
  stats(N, fmt(%12.0fc)) ///
  title("Weather Statistics during the Growing Season and Historical Comparison")
eststo clear


********************************************************************************
*** Table: Temperature bins box plots 
********************************************************************************

preserve 
	keep sh525less sh525 sh530 sh535plus sh525less_zscore sh525_zscore sh530_zscore sh535plus_zscore
	rename sh525less sh1 
	rename sh525 sh2
	rename sh530 sh3 
	rename sh535plus sh4 

	rename sh525less_zscore shz1 
	rename sh525_zscore shz2
	rename sh530_zscore shz3 
	rename sh535plus_zscore shz4 

	tempfile sh
	save sh, replace 
	forvalues j=1/4 {
		use sh
		keep sh`j' shz`j'
		rename sh`j' sh
		rename shz`j' shz
		gen bin = `j'
		tempfile sh`j'
		save sh`j', replace 
	}
	use sh1, clear
	append using sh2 
	append using sh3 
	append using sh4 
	label define bins 1 "<25C" 2 "25-30C" 3 "30-35C" 4 ">35C"
	label values bin bins 
	label var bin "Percentage of growing season hours"
	label var sh "Percentage of growing season hours"
	label var shz "Z-score"
	graph box sh, over(bin, label(labsize(large))) ylabel(,labsize(large)) ///
	ytitle(,size(large)) box(1,color(green) lcolor(black)) ///
	marker(1, mcolor(green)) yline(0) saving(bins, replace) 
	graph box shz, over(bin,  label(labsize(large))) ylabel(,labsize(large)) ///
	ytitle(,size(large)) box(1,color(green) lcolor(black)) ///
	marker(1, mcolor(green)) yline(0) saving(binszscore, replace) xsize(medium)
	graph save binszscore.gph, replace 
	graph close 
	graph combine bins.gph binszscore.gph
	graph export "${output}\Figures\fig4_temperature_bins_distribution.pdf", replace
restore 

********************************************************************************
*** Table: Sample selection (still need to label vars)
********************************************************************************

preserve 
	use "${data}/final_dataset.dta", clear 
	global dstats_sample "d_pyfull d_pyfullnon0"
	global dstats_sample "$dstats_sample d_plantedareafull d_plantedareafullnon0"
	global dstats_sample "$dstats_sample d_moduleany d_fullnon0_aux d_final"
	
	label var d_pyfull "Observations with this module (\% of full sample) "
	label var d_pyfullnon0 "Observations with module and strictly positive output (\% of full sample)"
	
	label var d_plantedareafull "Observations with this module (\% of full sample)"
	label var d_plantedareafullnon0 "Observations with module and strictly positive planted area (\% of full sample)"
	
	label var d_moduleany "Observations with at least one of the two modules (\% of full sample)"
	label var d_fullnon0_aux "Observations with strictly positive output and planted area (\% of full sample)"
	label var d_final "Final sample (% of full sample)"

	estpost summ $dstats_sample 
	est store dstats_all
	estpost summ $dstats_sample if year == 2010 
	est store dstats2010
	estpost summ $dstats_sample if year == 2012 
	est store dstats2012
	estpost summ $dstats_sample if year == 2015
	est store dstats2015

	esttab dstats_all dstats2010 dstats2012 dstats2015 , ///
		main(mean) /// aux(sd)
		refcat(	d_pyfull "Agrciultural Production Module" ///
				d_plantedareafull "Planted Area Module" ///
				d_moduleany " " ///
				, nolabel) ///
		mtitle("All" "2010" "2012" "2015" "2018") label ///
		stats(N, fmt(%12.0fc)) nonotes compress nogap lines 
	eststo clear 
	
	esttab dstats_all dstats2010 dstats2012 dstats2015 using "${output}/Tables/tab3_sample_selections_stats.tex", replace ///
		main(mean) /// aux(sd)
		refcat(	d_pyfull "Agrciultural Production Module" ///
				d_plantedareafull "Planted Area Module" ///
				d_moduleany " " ///
				, nolabel) ///
		mtitle("All" "2010" "2012" "2015" "2018") label ///
		stats(N, fmt(%12.0fc)) nonotes compress nogap lines ///
		title("Sample Selection")
	eststo clear 
restore 

********************************************************************************
*** Figure: Mean growing season temperatures in Nigeria over the period 1983-2015 
********************************************************************************

preserve
	import delimited "${root_dir}/Data/temp_bins_hours.csv", clear 
	gen temp_avg = 0 
	forvalues t = 15/49 {
		replace temp_avg = temp_avg + h`t'*`t'
	}
	replace temp_avg = temp_avg/(182*24)
	summ temp_avg 
	collapse (mean) temp_avg, by(year)
	label var temp_avg "Average temperature (C)"
	label var year "Year"
	graph twoway (line temp_avg year, color(green) lwidth(thick)), /// 
	xscale(range(1983 2014)) xlabel(1985(5)2015, labsize(medium)) ylabel(, labsize(medium))  ///
	ytitle(,size(medium)) xtitle(,size(medium))
	graph export "${output}/Figures/fig5_mean_temperature.pdf", replace
restore 

********************************************************************************
*** Table 4: Effect of Temperatures on Agricultural Productivity (to do: Oster's delta)
*******************************************************************************

eststo t41: reg2hdfespatial ln_output_2010usd $tempbins_linear $ltempbins_linear $precip $inputs_ln $inputs_dummies $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(100) lag(10)
qui estadd local farm Yes
qui estadd local year Yes
qui estadd local hh Yes
qui estadd local inputs Yes
qui estadd local interviewmonth Yes 
estimates store t41

eststo t42: reg2hdfespatial ln_yield_per_ha $tempbins_linear $ltempbins_linear $precip $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(100) lag(10)
qui estadd local farm Yes
qui estadd local year Yes
qui estadd local hh Yes
qui estadd local inputs No 
qui estadd local interviewmonth Yes 
estimates store t42

eststo t43: reg2hdfespatial ln_output_2010usd $tempbins_linear $ltempbins_linear $precip $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(100) lag(10)
qui estadd local farm Yes
qui estadd local year Yes
qui estadd local hh Yes
qui estadd local inputs No 
qui estadd local interviewmonth Yes 
estimates store t43

eststo t44: reg2hdfespatial ln_ha_planted $tempbins_linear $ltempbins_linear $precip $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(100) lag(10)
qui estadd local farm Yes
qui estadd local year Yes
qui estadd local hh Yes
qui estadd local inputs No
qui estadd local interviewmonth Yes 
estimates store t44

esttab t41 t42 t43 t44 ///
 ,  keep(sh535plus lsh535plus) ///
 order(sh535plus lsh535plus) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("TFP" "ln(yield)" "ln(value of output)" "ln(planted area)") ///
stats(farm year interviewmonth hh inputs N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Input Controls" "Observations" "R-squared" ))  ///
title("Other Outcomes") ///
star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .")

esttab t41 t42 t43 t44 using "${output}/Tables/tab4_main_results.tex" ///
 , replace keep(sh535plus lsh535plus) ///
 order(sh535plus lsh535plus) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("TFP" "ln(yield)" "ln(value of output)" "ln(planted area)") ///
stats(farm year interviewmonth hh inputs N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Input Controls" "Observations" "R-squared" ))  ///
title("Effect of Temperatures on Agricultural Productivity") ///
star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .")

********************************************************************************
*** Table 5: Adaptation to high temperatures: Extensive margin of input use (done) 
********************************************************************************  

*** Extensive margin 

local j = 0 
foreach var in labor_family labor_hired inorgfert_kg urea_kg npk_kg herb_kg pest_kg  {
	local j = `j' + 1 
	eststo t5`j': reg2hdfespatial `var'_ext $tempbins_linear $ltempbins_linear $precip $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(100) lag(10)
	qui estadd local farm Yes
	qui estadd local year Yes
	qui estadd local hh Yes
	qui estadd local interviewmonth Yes 
	estimates store t5`j'
}

esttab t51 t52 t53 t54 t55 t56 t57, ///
  keep(sh535plus lsh535plus) ///
  order(sh535plus lsh535plus) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("Family labor" "Hired labor" "Inorganic fertilizer" "Urea" "NPK" "Herbicide" "Pesticide") ///
stats(farm year interviewmonth hh N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" )) ///
  title("Extensive margin of input use") ///
  star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .") 

esttab t51 t52 t53 t54 t55 t56 t57 using "${output}/Tables/tab5_input_use_response.tex", replace ///
  keep(sh535plus lsh535plus) ///
  order(sh535plus lsh535plus) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("Family labor" "Hired labor" "Inorganic fertilizer" "Urea" "NPK" "Herbicide" "Pesticide") ///
stats(farm year interviewmonth hh N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" )) ///
  title("Adaptation to High Temperatures: Extensive Margin of Input Use") ///
  star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .") 
 
********************************************************************************
**** Table 6: Adaptation to High Temperatures: Land Allocation Strategies
********************************************************************************  

global planted_mix "ha_planted_mono_pc ha_planted_mixed_pc"
local j=0
foreach var in $planted_mix {
	local j=`j'+1
	eststo t6a`j': reg2hdfespatial `var' $tempbins_linear $ltempbins_linear $precip $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(100) lag(10)
	qui estadd local farm Yes
	qui estadd local year Yes 
	qui estadd local hh Yes
	qui estadd local interviewmonth Yes 
	summ `var'  if ha_planted_mono_pc+ha_planted_mixed_pc>=0.999 & ha_planted_mono_pc+ha_planted_mixed_pc<=1.001
	local meanoutcomeaux = round(`r(mean)',0.01)
	estadd local meanoutcome `meanoutcomeaux' 
}

forvalues j=1/8 {
	eststo t6b`j': reg2hdfespatial sharepymc`j' $tempbins_linear $ltempbins_linear $precip $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(100) lag(10)
	qui estadd local farm Yes
	qui estadd local year Yes 
	qui estadd local hh Yes
	qui estadd local interviewmonth Yes 
	summ sharepymc`j'
	local meanoutcomeaux = round(`r(mean)',0.01)
	estadd local meanoutcome `meanoutcomeaux' 	
}

esttab t6a2 t6b1 t6b2 t6b3 t6b4 t6b5 t6b6 t6b7 t6b8,  ///
 keep(sh535plus lsh535plus) /// 
 order(sh535plus lsh535plus) se substitute(\_ _) label ///
 mtitle("Mixed cropping" "Cassava" "Millet" "Sorghum" "Cowpeas" "Maize" "Rice" "Ground nuts" "Soybean") ///
stats(meanoutcome farm hh interviewmonth N r2, labels("Mean outcome" "Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" ))  ///
star(* 0.10 ** 0.05 *** 0.01)    addnote("{\it Notes:} .")  title("Table N: Ratio of total value of agricultural prodution")

esttab t6a2 t6b1 t6b2 t6b3 t6b4 t6b5 t6b6 t6b7 t6b8 using "${output}/Tables/tab6_land_allocation_strategies.tex", replace  ///
 keep(sh535plus lsh535plus) /// 
 order(sh535plus lsh535plus) se substitute(\_ _) label ///
 mtitle("Mixed cropping" "Cassava" "Millet" "Sorghum" "Cowpeas" "Maize" "Rice" "Ground nuts" "Soybean") ///
stats(meanoutcome farm hh interviewmonth N r2, labels("Mean outcome" "Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" ))  ///
star(* 0.10 ** 0.05 *** 0.01)    addnote("{\it Notes:} .")  title(" Adaptation to High Temperatures: Land Allocation Strategies")


********************************************************************************
*** Figure 6: Falsification Tests 
********************************************************************************
ssc install shufflevar
preserve 
	global tempbins_linear_shuffle sh525less_shuffle sh530_shuffle sh535plus_shuffle lsh525less_shuffle
	global precip_shuffle rf_shuffle rf2_shuffle rf_lag_shuffle rf2_lag_shuffle
	gen coeff_scramble_hhid = . 

	gen ln_yield = ln_yield_per_ha
	gen ln_area = ln_ha_planted

	set seed 123456789
	foreach var in yield area {
		foreach level in year {
			gen coeff_scramble_`var'_`level' = .
			forvalues j=1/1000 {
				shufflevar $tempbins_linear $ltempbins_linear $precip , joint cluster(`level')
				qui reg ln_`var' $tempbins_linear_shuffle  $precip_shuffle $hhchar $soilhalf interviewmonth1-interviewmonth9
				display _b[sh535plus_shuffle]
				matrix b=e(b)
				display b[1,3]
				replace coeff_scramble_`var'_`level' = b[1,3] in `j'
			}
		}
	}
	
	histogram coeff_scramble_yield_year, fraction xscale(range(-0.09 0.03)) ///
	color(green) lcolor(black) xlabel(-0.09(0.01)0.03) xline(-0.077) xtitle("Coefficient") 
	graph export "${output}/Figures/fig6_coeff_scramble_yield_year.pdf", replace
	graph save "${output}/Figures/fig6_coeff_scramble_yield_year.gph", replace
	histogram coeff_scramble_area_year, fraction xlabel(-0.02(0.01)0.05) ///
	color(green) lcolor(black) xline(0.0429) xtitle("Coefficient")
	graph export "${output}/Figures/fig6_coeff_scramble_area_year.pdf", replace
	graph save "${output}/Figures/fig6_coeff_scramble_area_year.gph", replace
restore 

********************************************************************************
*** Table: Spatial-HAC Standard Errors - Sensitivity
********************************************************************************

foreach temp_cutoff in 0 5 10 {
	foreach dist_cutoff in 0 100 250 {
		display "$tempbins_linear"

		eststo t71: reg2hdfespatial ln_output_2010usd $tempbins_linear $ltempbins_linear $precip $inputs_ln $inputs_dummies $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(`dist_cutoff') lag(`temp_cutoff')
		qui estadd local farm Yes
		qui estadd local year Yes
		qui estadd local hh Yes
		qui estadd local inputs Yes
		qui estadd local interviewmonth Yes 
		estimates store t32_b

		eststo t72: reg2hdfespatial ln_yield_per_ha $tempbins_linear $ltempbins_linear $precip $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(`dist_cutoff') lag(`temp_cutoff')
		qui estadd local farm Yes
		qui estadd local year Yes
		qui estadd local hh Yes
		qui estadd local inputs No 
		qui estadd local interviewmonth Yes 
		estimates store t31

		eststo t73: reg2hdfespatial ln_output_2010usd $tempbins_linear $ltempbins_linear $precip $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(`dist_cutoff') lag(`temp_cutoff')
		qui estadd local farm Yes
		qui estadd local year Yes
		qui estadd local hh Yes
		qui estadd local inputs No 
		qui estadd local interviewmonth Yes 
		estimates store t32
		
		eststo t74: reg2hdfespatial ln_ha_planted $tempbins_linear $ltempbins_linear $precip $hhchar $soilhalf interviewmonth1-interviewmonth9, timevar(year) panelvar(hhid) lat(lat_dd_mod) lon(lon_dd_mod) distcutoff(`dist_cutoff') lag(`temp_cutoff')
		qui estadd local farm Yes
		qui estadd local year Yes
		qui estadd local hh Yes
		qui estadd local inputs No
		qui estadd local interviewmonth Yes 
		estimates store t33

esttab t71 t72 t73 t74 ///
 ,  keep($tempbins_linear $ltempbins_linear) ///
 order($tempbins_linear $ltempbins_linear) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("ln(yield)" "ln(value of output)" "TFP" "ln(planted area)") ///
stats(farm year interviewmonth hh inputs N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Input Controls" "Observations" "R-squared" ))  ///
title("Other Outcomes") ///
star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .")

esttab t71 t72 t73 t74 using "${output}/Tables/tab7_sensitivity_`temp_cutoff'_`dist_cutoff'.tex", replace ///
  keep($tempbins_linear $ltempbins_linear) ///
  order($tempbins_linear $ltempbins_linear) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("ln(yield)" "ln(value of output)" "TFP" "ln(planted area)") ///
stats(farm year interviewmonth hh inputs N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Input Controls" "Observations" "R-squared" )) ///
  title("Other Outcomes - Conley SEs with temporal lag of `temp_cutoff' years and distance cut-off of `dist_cutoff' km") ///
  star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .") 
	}
}


********************************************************************************
*** Table: Predicted Effects of a Uniform Warming Scenario
********************************************************************************

* Generate counterfactual hours per 1-degree-Celsius temperature bins  
foreach hours in h lh {
	gen `hours'15_uws=0
	forvalues j = 15/49 {
		local k = `j'+1
		capture drop  `hours'`k'_uws
		gen `hours'`k'_uws= `hours'`j'
	}
}

* Generate counterfactual temperature bins 
foreach share in sh lsh {
	forvalues temp=15(5)45 {
		local temp4 = `temp' + 4
		capture drop `share'5`temp'_uws
		if "`share'"=="sh" {
			egen `share'5`temp'_uws=rowtotal(h`temp'_uws-h`temp4'_uws)
		}
		if "`share'"=="lsh" {
			egen `share'5`temp'_uws=rowtotal(lh`temp'_uws-lh`temp4'_uws)
		}
		replace `share'5`temp'_uws=100*`share'5`temp'_uws/(182*24)
		summ `share'5`temp'_uws,d 
		}
	gen `share'540plus_uws = `share'540_uws + `share'545_uws
	gen `share'535plus_uws = `share'535_uws + `share'540plus_uws
	gen `share'525less_uws = `share'515_uws + `share'520_uws
	summ `share'535plus_uws, d 
	summ `share'530_uws, d 
	summ `share'525_uws, d 
	summ `share'525less_uws, d 
}

* Generate variables with change in temperature bins from baseline to counterfactual 
foreach var in $tempbins_linear $ltempbins_linear {
	display "`var'"
	capture drop `var'_delta
	gen `var'_delta = `var'_uws-`var'
}

* Call regression results 
forvalues j=1/4 {
	estimates restore t4`j'
}
esttab t41 t42 t43 t44 ///
	,  keep($tempbins_linear) label ///
	mtitle("ln(yield)" "ln(value of output)" "TFP" "ln(planted area)") ///
	stats(farm year interviewmonth hh inputs N r2, ///
	labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Input Controls" "Observations" "R-squared" ))  ///
	title("Main Outcomes") ///
	star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .")

* Get counterfactual responses 
foreach j in 1 2 3 4 {
	estimates restore t4`j' 
	matrix coeff = e(b)
	local k = 0 
	foreach var in $tempbins_linear $ltempbins_linear {
		local k = `k'+1
		local coeff_`j'_`k'=_b["`var'"]
		display "`coeff_`j'_`k''"
	}
	capture drop change_`j'
	gen change_`j' = `coeff_`j'_1'*sh525less_delta
	local m = 1
	foreach var in sh530 sh535plus $ltempbins_linear {
		local m = `m'+1 
		replace change_`j' = change_`j'+`coeff_`j'_`m''*`var'_delta
	}
}

* Show results 
foreach j in 1 2 3 4 {
	summ  change_`j' [w=ha_planted]
}

* Drop created variables 
drop h15_uws-change_4 

********************************************************************************
*** Appendix Table: Contemporaneous temperature and sample selection
********************************************************************************

global hhcharacteristics "age age2 edu female hhsize_num"
global tempbins "sh525less sh530 sh535plus"
global rf "rf rf2"
global tempbins_mean "sh525less_mean sh530_mean sh535plus_mean"
global rf_mean "rf_mean rf2_mean"
global tempbins_zscore "sh525less_zscore sh530_zscore sh535plus_zscore"
global rf_zscore "rf_zscore rf2_zscore"

preserve 
	use "${data}/final_dataset.dta", clear 
	eststo ta11: probit d_final  $tempbins $rf $hhcharacteristics i.ssa_aez09 if year==2010, vce(cluster lga)
	qui estadd local year 2010
	eststo ta12: probit d_final  $tempbins $rf  $hhcharacteristics i.ssa_aez09 if year==2012, vce(cluster lga)
	qui estadd local year 2012
	eststo ta13: probit d_final  $tempbins $rf  $hhcharacteristics i.ssa_aez09 if year==2015, vce(cluster lga)
	qui estadd local year 2015
	
	esttab ta11 ta12 ta13 ///
	 ,  keep($tempbins $rf  $hhcharacteristics) ///
	 order($tempbins $rf  $hhcharacteristics)  se label ///
	stats(year N r2_p , labels("Year" "Observations" "R-squared")) collabels("2010" "2012" "2015") ///
	star(* 0.10 ** 0.05 *** 0.01)   nonotes
	
	esttab ta11 ta12 ta13 using "${output}/Tables/taba1_sample_selection_current_temp.tex", replace ///
	keep($tempbins $rf  $hhcharacteristics) ///
	 order($tempbins $rf  $hhcharacteristics)  se label ///
	stats(year N r2_p , labels("Year" "Observations" "R-squared")) collabels("2010" "2012" "2015") ///
	star(* 0.10 ** 0.05 *** 0.01)   nonotes ///
	title("Contemporaneous temperature and sample selection")
restore 

********************************************************************************
*** Appendix Table: Long-term temperature and sample selection
********************************************************************************

preserve 
	use "${data}/final_dataset.dta", clear 
	eststo ta21: probit d_final   $tempbins_mean $rf_mean $tempbins_zscore $rf_zscore $hhcharacteristics i.ssa_aez09 if year==2010, vce(cluster lga)
	qui estadd local year 2010
	eststo ta22: probit d_final   $tempbins_mean $rf_mean $tempbins_zscore $rf_zscore $hhcharacteristics i.ssa_aez09 if year==2012, vce(cluster lga)
	qui estadd local year 2012
	eststo ta23: probit d_final   $tempbins_mean $rf_mean $tempbins_zscore $rf_zscore $hhcharacteristics i.ssa_aez09 if year==2015, vce(cluster lga)
	qui estadd local year 2015
	
	esttab ta21 ta22 ta23 ///
	 ,  keep( $tempbins_mean $rf_mean $tempbins_zscore $rf_zscore $hhcharacteristics) ///
	 order( $tempbins_mean $rf_mean $tempbins_zscore $rf_zscore $hhcharacteristics)  se label ///
	stats(year N r2_p , labels("Year" "Observations" "R-squared")) collabels("2010" "2012" "2015") ///
	star(* 0.10 ** 0.05 *** 0.01)   addnote("{\i Notes:} .")

	esttab ta21 ta22 ta23 using "${output}/Tables/taba2_sample_selection_lagged_temp.tex", replace ///
	 keep( $tempbins_mean $rf_mean $tempbins_zscore $rf_zscore $hhcharacteristics) ///
	 order( $tempbins_mean $rf_mean $tempbins_zscore $rf_zscore $hhcharacteristics)  se label ///
	stats(year N r2_p , labels("Year" "Observations" "R-squared")) collabels("2010" "2012" "2015") ///
	star(* 0.10 ** 0.05 *** 0.01)   addnote("{\i Notes:} .")
restore 

********************************************************************************
*** Appendix Table A3: Effect of Temperatures on Agricultural Productivity - Full Results
********************************************************************************

esttab t41 t42 t43 t44 ///
 ,  keep($tempbins_linear $ltempbins_linear $precip $hh_chars) ///
 order($tempbins_linear $ltempbins_linear $precip $hh_chars) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("TFP" "ln(yield)" "ln(value of output)" "ln(planted area)") ///
stats(farm year interviewmonth hh inputs N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Input Controls" "Observations" "R-squared" ))  ///
title("Other Outcomes") ///
star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .")

esttab t41 t42 t43 t44 using "${output}/Tables/taba3_main_results_full.tex", replace ///
 keep($tempbins_linear $ltempbins_linear $precip $hhchar) ///
 order($tempbins_linear $ltempbins_linear $precip $hhchar) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("TFP" "ln(yield)" "ln(value of output)" "ln(planted area)") ///
stats(farm year interviewmonth hh inputs N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Input Controls" "Observations" "R-squared" ))  ///
title("Effect of Temperatures on Agricultural Productivity - Full Results") ///
star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .")

********************************************************************************
*** Appendix Table A4: Adaptation to High Temperatures: Extensive Margin of Input Use - Full Results
********************************************************************************

esttab t51 t52 t53 t54 t55 t56 t57, ///
  keep($tempbins_linear $ltempbins_linear $precip $hhchar) ///
  order($tempbins_linear $ltempbins_linear $precip $hhchar) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("Family labor" "Hired labor" "Inorganic fertilizer" "Urea" "NPK" "Herbicide" "Pesticide") ///
stats(farm year interviewmonth hh N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" )) ///
  title("Adaptation to High Temperatures: Extensive Margin of Input Use - Full Results") ///
  star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .") 

esttab t51 t52 t53 t54 t55 t56 t57 using "${output}/Tables/taba4_input_use_full.tex", replace ///
  keep($tempbins_linear $ltempbins_linear $precip $hhchar) ///
  order($tempbins_linear $ltempbins_linear $precip $hhchar) se substitute("<" "$<$" ">" "$>$") label ///
 mtitle("Family labor" "Hired labor" "Inorganic fertilizer" "Urea" "NPK" "Herbicide" "Pesticide") ///
stats(farm year interviewmonth hh N r2, labels("Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" )) ///
  title("Adaptation to High Temperatures: Extensive Margin of Input Use - Full Results") ///
  star(* 0.10 ** 0.05 *** 0.01)   addnote("{\it Notes:} .") 
 
  
********************************************************************************
*** Appendix Table A5: Adaptation to High Temperatures: Land Allocation Strategies - Full Results
********************************************************************************

esttab t6a2 t6b1 t6b2 t6b3 t6b4 t6b5 t6b6 t6b7 t6b8,  ///
 keep($tempbins_linear $ltempbins_linear $precip $hhchar) /// 
 order($tempbins_linear $ltempbins_linear $precip $hhchar) se substitute(\_ _) label ///
 mtitle("Cassava" "Millet" "Sorghum" "Cowpeas" "Maize" "Rice" "Ground nuts" "Soybean") ///
stats(meanoutcome farm year hh interviewmonth N r2, labels("Mean outcome" "Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" ))  ///
star(* 0.10 ** 0.05 *** 0.01)    addnote("{\it Notes:} .")  title("Adaptation to High Temperatures: Land Allocation Strategies - Full Results")

esttab t6a2 t6b1 t6b2 t6b3 t6b4 t6b5 t6b6 t6b7 t6b8 using "${output}/Tables/taba5_land_allocation_full.tex", replace ///
 keep($tempbins_linear $ltempbins_linear $precip $hhchar) /// 
 order($tempbins_linear $ltempbins_linear $precip $hhchar) se substitute(\_ _) label ///
 mtitle("Cassava" "Millet" "Sorghum" "Cowpeas" "Maize" "Rice" "Ground nuts" "Soybean") ///
stats(meanoutcome farm year hh interviewmonth N r2, labels("Mean outcome" "Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" ))  ///
star(* 0.10 ** 0.05 *** 0.01)    addnote("{\it Notes:} .")  title("Adaptation to High Temperatures: Land Allocation Strategies - Full Results")


********************************************************************************
*** Appendix Table A6: Adaptation to High Temperatures: Land Allocation Strategies - Full Results
********************************************************************************

esttab t6a1 t6a2,  ///
 keep($tempbins_linear $ltempbins_linear $precip $hhchar) /// 
 order($tempbins_linear $ltempbins_linear $precip $hhchar) se substitute(\_ _) label ///
 mtitle("Mono cropping" "Mixed cropping") ///
stats(meanoutcome farm year hh interviewmonth N r2, labels("Mean outcome" "Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" ))  ///
star(* 0.10 ** 0.05 *** 0.01)    addnote("{\it Notes:} .")  title("Adaptation to High Temperatures: Land Allocation Strategies - Full Results")

esttab t6a1 t6a2 using "${output}/Tables/taba6_land_allocation_full.tex", replace ///
 keep($tempbins_linear $ltempbins_linear $precip $hhchar) /// 
 order($tempbins_linear $ltempbins_linear $precip $hhchar) se substitute(\_ _) label ///
 mtitle("Mono cropping" "Mixed cropping") ///
stats(meanoutcome farm year hh interviewmonth N r2, labels("Mean outcome" "Farm FE" "Year FE" "Interview Month FE" "HH Characteristics" "Observations" "R-squared" ))  ///
star(* 0.10 ** 0.05 *** 0.01)    addnote("{\it Notes:} .")  title("Adaptation to High Temperatures: Land Allocation Strategies - Full Results")

